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(1)  Objectives: 

Recently,  the  research  on  the  flapping  micro  air  vehicle  which  imitates  living  creatures  is  being  actively 
carried  out.  The  objective  of  this  research  is  the  definition  of  the  work  flow  for  the  design  of  the  flapping 
wing.  Generally,  the  best  performance  and  the  minimum  mass  should  be  pursued  in  the  design  of  the 
flapping  wing  micro  air  vehicle.  The  flapping  wing  structure  should  have  some  rigidity  which  is  kept  with 
minimum  mass.  Therefore,  path  optimization  is  performed  for  maximizing  the  performance  and  dynamic 
topology  optimization  is  performed  for  determining  the  distribution  of  the  reinforcement.  Reinforcement 
within  the  wing  is  determined  to  have  certain  rigidity  with  minimum  mass.  Dynamic  topology  optimization 
is  performed  based  on  the  optimum  values  of  path  optimization. 

(2)  Status  of  effort: 

This  research  starts  from  a  well  defined  flapping  wing  which  is  proposed  by  Jones  and  Platzer  in  2006. 
Path  optimization  is  performed  for  a  flapping  wing.  The  objective  function  to  be  maximized  is  a  multi¬ 
objective  function  which  is  composed  of  the  thrust  coefficient  and  the  propulsive  efficiency.  In  order  to 
evaluate  the  objective  function,  computational  fluid  dynamics  (CFD)  analysis  is  used.  A  gradient  based 
optimization  is  extremely  expensive  because  of  CFD  analysis.  Therefore,  a  surrogate  model  is  utilized  and 
the  Kriging  method  is  employed  to  establish  the  surrogate  model.  The  design  points  (samples)  are  defined  by 
an  orthogonal  array.  A  genetic  algorithm  is  selected  to  solve  the  multi-objective  function  problem.  For 
structural  design,  3D  unsteady  CFD  analysis  is  performed  with  the  optimum  values  from  the  path 
optimization.  The  3D  unsteady  CFD  analysis  results  supply  the  pressure  distribution  in  the  time  domain  and 
this  dynamic  pressure  distribution  is  used  as  the  external  loads  for  dynamic  analysis  of  the  flapping  wing 
structure.  Dynamic  topology  optimization  is  performed  to  find  the  distribution  of  reinforcement.  A  dynamic 
topology  optimization  method  is  developed  by  modification  of  the  Equivalent  Static  Loads  method  for  non 
linear  static  response  Structural  Optimization  (ESLSO). 

(3)  Abstract: 

The  flapping  wing  of  a  micro  air  vehicle  is  optimized  to  enhance  performance  while  some  rigidity  is  kept 
with  minimum  mass.  A  work  flow  for  the  design  of  the  flapping  wing  is  defined.  The  performances  to  be 
enhanced  are  thrust  coefficient  and  propulsive  efficiency.  The  flapping  kinematics  of  the  flapping  wing  is 
determined  by  solving  a  path  optimization  problem  which  maximizes  the  performances.  The  optimization 


process  is  carried  out  based  on  a  well  defined  surrogate  model.  The  surrogate  model  is  made  from  the  results 
of  two-dimensional  fluid  dynamic  analysis.  The  Kriging  method  is  employed  to  establish  the  surrogate 
model  and  a  genetic  algorithm  is  utilized  for  the  multi-objective  function  problem.  Dynamic  topology 
optimization  is  performed  to  find  the  distribution  of  reinforcement.  Certain  rigidity  can  be  kept  by  the  results 
of  topology  optimization.  A  dynamic  topology  optimization  method  is  developed  by  modification  of  the 
equivalent  static  loads  method  for  non  linear  static  response  structural  optimization.  Three-dimensional 
computational  fluid  dynamic  analysis  is  performed  based  on  the  optimum  values  of  the  path  optimization  to 
evaluate  the  external  loads  for  the  topology  optimization  process.  It  is  found  that  the  topology  results  are 
quite  similar  to  the  practical  product.  The  process  of  the  defined  work  flow  is  materialized  by  interfacing 
various  software  systems. 
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Abstract 


The  flapping  wing  of  a  micro  air  vehicle  is  optimized  to  enhance  performance  while  some  rigidity  is  kept 
with  minimum  mass.  A  work  flow  for  the  design  of  the  flapping  wing  is  defined.  The  performances  to  be 
enhanced  are  thrust  coefficient  and  propuisive  efficiency.  The  flapping  kinematics  of  the  flapping  wing  is 
determined  by  solving  a  path  optimization  problem  which  maximizes  the  performances.  The  optimization 
process  is  carried  out  based  on  a  well  defined  surrogate  model.  The  surrogate  model  is  made  from  the  results 
of  two-dimensional  fluid  dynamic  analysis.  The  Kriging  method  is  employed  to  establish  the  surrogate  model 
and  a  genetic  algorithm  is  utilized  for  the  multi-objective  function  problem.  Dynamic  topology  optimization 
is  performed  to  find  the  distribution  of  reinforcement.  Certain  rigidity  can  be  kept  by  the  results  of  topology 
optimization.  A  dynamic  topology  optimization  method  is  developed  by  modification  of  the  equivalent  static 
loads  method  for  non  linear  static  response  structural  optimization.  Three-dimensional  computational  fluid 
dynamic  analysis  is  performed  based  on  the  optimum  values  of  the  path  optimization  to  evaluate  the  external 
loads  for  the  topology  optimization  process.  It  is  found  that  the  topology  results  are  quite  similar  to  the 
practical  product.  The  process  of  the  defined  work  flow  is  materialized  by  interfacing  various  software 
systems. 


Nomenclature 

AR  =  aspect  ratio 

a  =  dimensionless  length  from  the  leading  edge  to  the  pivot  point 


b  =  design  variable  vector 

C(b)  =  viscous  damping  matrix 

CFD  =  computational  fluid  dynamics 

C/j  =  drag  coefficient 

Ci  =  lift  coefficient 

Ct  =  thrust  coefficient 

CM  =  pitch  moment  coefficient 

CD  =  time  average  drag  coefficient 

C  =  time  average  lift  coefficient 

C  =  time  average  thrust  coefficient 

c  =  chord,  m 

E  =  Young’s  modulus 

Fx  =  x  component  of  the  resulting  aerodynamics  force  acting  on  the  airfoil,  N 

Fy  =  y  component  of  the  resulting  aerodynamics  force  acting  on  the  airfoil,  N 

f  =  frequency.  Hertz 

f  (s)  =  equivalent  static  loads  vector 

h  =  reduced  plunging  amplitude  with  respect  to  the  chord 

k  =  reduced  frequency 

K(b)  =  stiffness  matrix 

M(b)  =  mass  matrix 

Ma  =  mach  number 

Re  =  Reynolds  number 

S  =  reference  area 

St  =  Strouhal  number 

T  =  period,  second 

u  =  far  field  flow  velocity,  m/s 
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total  volume  of  the  structure 
dynamic  displacement  vector 
static  displacement  vector 

phase  angle  of  pitching  motion  leading  plunging  motion 
density  of  fluid  around  the  airfoil,  kg/m3 
density  of  the  material,  kg/m  ' 

angle  between  the  chord  and  the  far  field  flow  speed  direction 

propulsive  efficiency 

angular  frequency,  rad/s 

specified  value  of  density 

percent  of  the  total  design  variable 

density  value  of  design  update 


I.  Introduction 

The  flapping  wing  MAV  (FWMAV),  inspired  by  birds,  bats,  insects,  fishes  and  whales,  has  been  getting  more  and 
more  attention  from  military  and  civilian  application  domains  since  the  micro  air  vehicle  (MAV)  was  generally 
defined  by  the  Defense  Advanced  Research  Projects  Agency  (DARPA)  in  1997  [1],  Some  unique  advantages  of 
FWMAVs,  such  as  high  maneuverability  and  hovering  capability,  make  the  research  in  this  area  quite  active.  Early 
studies,  such  as  the  Garrick  theory  [2],  focused  on  simplified  models.  As  research  evolved,  practical  and 
sophisticated  models  are  employed.  The  continuous  development  of  computational  methods  makes  it  possible  to 
simulate  a  real  insect  flying  in  the  air  [3]  and  some  FWMAVs  were  fabricated  during  the  past  years.  For  example, 
Jones  and  Platzer  proposed  an  unconventional  biplane  flapping  wing  MAV  in  2006  [4].  In  this  model,  the  thrust  is 
generated  by  the  biplane  pair  of  the  two  trailing  flapping  wings  and  the  lift  is  generated  by  the  front  stationary  wing 
as  shown  in  Fig.  1.  The  present  work  is  closely  related  to  this  design. 

From  the  perspective  of  engineering,  the  best  performance  and  the  minimum  mass  should  be  pursued  in  the  design 
of  the  FWMAV.  The  flapping  path  should  be  determined  to  maximize  the  performance.  Recent  experimental  and 
computational  achievements  showed  that  the  flapping  performance,  such  as  the  thrust  coefficient  and  the  propulsive 


efficiency,  depends  on  the  flapping  kinematics  significantly.  Anderson  showed  that  oscillating  foil  could  have  a 
very  high  propulsive  efficiency,  as  high  as  87%,  under  specific  combinations  of  the  flapping  kinematics  by  water 
tunnel  experiments  [5].  Pesavento  and  Wang  found  that  optimized  flapping  wing  motions  could  save  up  to  27%  of 
the  aerodynamic  power  required  by  the  optimal  steady  flight  [6],  In  2005,  Tuncer  and  Kaya  employed  the  steepest 
ascent  search  algorithm  to  maximize  a  linear  combination  of  the  maximum  thrust  and  the  propulsive  efficiency  for 
flapping  airfoils  [7].  In  2007,  Kaya  and  Tuncer  performed  optimization  with  a  NURBS  (Non-Uniform  Rational  B- 
Splines,  NURBS)  flapping  path  instead  of  a  sinusoidal  flapping  path  in  their  previous  paper,  and  a  gradient  based 
search  algorithm  was  employed  to  find  the  optimum  [8].  In  2009,  Kaya,  Tuncer,  Jones  and  Platzer,  extended  their 
optimization  scheme  to  a  biplane  configuration  [9].  Meanwhile,  Soueid,  Guglielmini,  Airiau  and  Bottaro  optimized 
the  motion  of  a  flapping  airfoil  using  sensitivity  functions.  In  Soueid’s  work,  the  objective  was  a  combination  of 
several  concerned  parameters,  such  as  the  time  average  thrust  coefficient,  the  average  power  input  and  the  average 
angle  of  attack  [10]. 

In  this  research,  path  optimization  is  performed  based  on  the  paradigm  of  the  previous  researches  [7];  however, 
different  optimization  methods  are  utilized  to  reduce  the  effort  for  path  optimization.  A  surrogate  model  is  utilized 
in  the  optimization  process  and  the  Kriging  method  is  employed  to  make  the  surrogate  model  [11].  The  surrogate 
model  is  established  by  using  2-dimensional  computational  fluid  dynamics  (2D  CFD)  analysis  because  3D  CFD 
analysis  in  the  optimization  process  is  uncontrollably  expensive.  The  surrogate  model  has  the  advantage  over  a 
gradient  based  optimization  in  that  sensitivity  analysis  is  not  necessary.  It  is  well  known  that  sensitivity  information 
from  CFD  analysis  is  extremely  expensive.  In  other  words,  the  surrogate  model  is  used  to  reduce  the  computational 
cost.  Since  the  objective  function  is  a  multi-objective  function  defined  by  the  thrust  coefficient  and  the  propulsive 
efficiency,  a  genetic  algorithm  is  utilized  [12].  Moreover,  the  genetic  algorithm  can  find  a  global  solution  while  a 
gradient  based  optimization  can  find  a  local  optimum. 

Topology  optimization  is  carried  out  to  keep  certain  rigidity  of  the  flapping  wing  with  minimum  mass.  It  is  noted 
that  there  are  not  many  studies  on  the  structural  design  of  the  flapping  wing.  Some  researchers  have  conducted 
topology  optimization  of  the  fixed  wing  of  MAV  [13].  In  this  research  a  topology  optimization  method  is  newly 
developed  for  the  flapping  wing.  Topology  optimization  determines  the  distribution  of  the  reinforcement  on  the 
flapping  wing.  The  distribution  of  the  reinforcement  can  be  regarded  as  the  vein  distribution  of  the  wing  which 
mimics  the  wing  of  an  insect.  We  need  pressure  distribution  on  the  flapping  wing.  Since  2D  CFD  analysis  does  not 


provide  pressure  distribution  on  the  entire  body,  3D  CFD  analysis  is  performed  based  on  the  optimum  values  of  path 
optimization.  The  pressure  from  the  3D  CFD  analysis  is  imposed  on  the  flapping  wing  as  the  external  loads  in  the 
topology  optimization  process.  It  is  noted  that  the  pressures  are  imposed  in  the  time  domain;  therefore,  we  need 
dynamic  topology  optimization.  Generally,  topology  optimization  is  performed  in  a  static  sense.  A  dynamic 
topology  optimization  method  is  developed  based  on  modification  of  a  method  called  the  Equivalent  Static  Loads 
method  for  non  linear  static  response  Structural  Optimization  (ESLSO)  [14-16]. 

The  problem  is  solved  by  interfacing  various  software  systems.  The  commercial  computational  fluid  dynamics 
(CFD)  system,  FLUENT  [17],  is  utilized  for  unsteady  aerodynamics  analysis.  NASTRAN  [18]  is  used  for  dynamic 
analysis  and  GENESIS  [19]  is  used  for  topology  optimization.  The  optimization  process  uses  the  genetic  algorithm 
in  MATLAB  [20].  An  in-house  C++  program  [21]  is  coded  to  link  the  systems.  The  research  results  are 
summarized  and  the  future  direction  for  the  design  of  the  flapping  wing  is  proposed. 


Fig.  1  Flapping-wing  MAV  model. 

II.  Flow  of  the  Present  Research 

This  work  starts  from  a  well  defined  flapping  wing  illustrated  in  Fig.  1.  The  biplane  pair  of  the  trailing  wings  is 
flapping  and  generates  the  thrust  force.  Path  optimization  is  performed  for  a  wing  of  the  two  flapping  wings.  The 
overall  work  flow  is  presented  in  Fig.  2.  An  optimization  problem  to  determine  the  path  is  formulated.  The  motion 
of  the  flapping  wing  consists  of  the  plunging  motion  and  the  pitching  motion.  The  parameters  for  the  path  are  the 
plunging  amplitude,  pitching  amplitude  and  the  phase  angle  between  the  two  motions,  and  these  parameters  are  used 


as  design  variables  in  path  optimization.  The  objective  function  to  be  maximized  is  a  multi-objective  function 


which  is  composed  of  the  thrust  coefficient  and  the  propulsive  efficiency.  CFD  analysis  is  required  to  evaluate  the 
objective  function.  A  gradient  based  optimization  process  needs  sensitivity  information  with  respect  to  design 
variables  and  it  is  well  known  that  sensitivity  analysis  with  CFD  analysis  is  extremely  expensive  [10].  Therefore,  a 
surrogate  model  is  utilized  in  this  research  and  the  Kriging  method  is  employed  to  establish  the  surrogate  model  [11]. 
Because  3D  analysis  is  highly  costly,  2D  unsteady  CFD  analysis  is  conducted  for  the  establishment  of  the  surrogate 
model. 

Responses  such  as  the  thrust  coefficient  and  the  propulsive  efficiency  are  evaluated  at  various  design  points.  The 
design  points  (samples)  are  defined  by  an  orthogonal  array.  The  surrogate  model  from  the  Kriging  method  is  used 
in  the  optimization  problem.  A  genetic  algorithm  is  selected  to  solve  the  multi-objective  function  problem  [12]. 
The  genetic  algorithm  may  not  be  practical  if  we  directly  use  CFD  analysis  because  the  algorithm  requires  many 
function  calculations.  Since  the  surrogate  model  is  utilized,  the  function  calculation  in  the  optimization  process  is 
almost  negligible. 

For  structural  design,  the  pressures  on  the  flapping  wing  are  required  because  they  should  be  used  as  external 
loads  on  the  structure.  Only  3D  CFD  analysis  can  provide  the  pressure  distribution  around  the  flapping  wing. 
Therefore,  a  computational  model  in  three  dimensions  is  made  and  3D  unsteady  CFD  analysis  is  performed  with  the 
optimum  values  from  the  path  optimization.  The  flapping  motion  is  dynamic  (transient)  and  we  need  the  pressure 
distribution  in  the  time  domain.  The  3D  unsteady  CFD  analysis  results  supplies  the  pressure  distribution  in  the  time 
domain  and  this  dynamic  pressure  distribution  is  used  as  the  external  loads  for  dynamic  analysis  of  the  flapping 
wing  structure. 

The  flapping  wing  structure  should  have  some  flexibility  as  well  as  certain  rigidity  [22-23],  and  a  small  amount  of 
rigidity  guarantees  flexibility.  Topology  optimization  is  performed  to  satisfy  this  condition  with  minimum  mass. 
The  reinforcement  in  the  flapping  wing  can  be  determined  in  the  topology  optimization  process.  The  reinforcement 
can  be  regarded  as  the  vein  distribution  of  the  wing  if  the  flapping  wing  mimics  that  of  an  insect.  Topology 
optimization  is  generally  carried  out  with  static  loads  [24].  Flowever,  the  flapping  wing  operates  in  the  dynamic 
(transient)  environment.  Therefore,  the  current  topology  optimization  methods  can  be  hardly  used  for  the  design  of 
the  flapping  wing.  A  topology  optimization  method  is  developed  to  handle  the  dynamic  pressure  distribution  in  the 
time  domain.  The  new  method  is  a  modified  version  of  the  Equivalent  Static  Loads  method  for  non  linear  static 


response  Structural  Optimization  (ESLSO)  [25].  ESLSO  has  been  extensively  used  for  size  and  shape  optimizations 
and  it  is  modified  for  topology  optimization  in  this  research. 

The  logic  for  the  current  work  may  not  be  perfect  because  of  the  cost  and  the  limitation  of  the  technical  problems. 
For  example,  if  optimization  with  3D  aeroelasticity  analysis  is  used,  the  path  and  the  structure  can  be  optimized 
simultaneously.  However,  this  approach  does  not  seem  to  be  realistic  at  this  moment.  Moreover,  the  pair  of  the 
flapping  should  be  considered  together.  The  future  direction  of  this  research  is  proposed  in  a  later  section. 


III.  Path  Optimization 


This  section  starts  with  the  description  of  the  flapping  motion.  The  involved  governing  equations  and  a  well 
validated  mesh  for  simulating  the  unsteady  flow  field  around  an  airfoil  are  also  described.  After  the  path 
optimization  problem  is  defined,  the  surrogate  model  for  the  Kriging  method  is  carefully  constructed.  Then,  the 
optimum  flapping  kinematic  parameters  are  obtained  by  a  genetic  algorithm.  At  the  end  of  this  section,  3D  CFD 
analysis  is  performed  under  these  optimum  parameters  for  a  wing  with  a  given  aspect  ratio.  Meanwhile,  the 
derivations  between  flapping  performances  of  the  2D  airfoil  and  ones  of  the  3D  wing  are  demonstrated. 

A.  Flapping  Wing  System 

The  flapping  motion,  including  the  plunging  and  pitching  motions,  can  be  described  as  follows: 

y(t)  =  he  cos(a>t)  (1) 


9(t)  =  90  cos(ot  +  (/>)  (2) 

where  y(t)  stands  for  the  plunging  motion,  9(t)  the  pitching  motion,  c  is  the  chord  length,  h  is  the  dimensionless 
plunging  amplitude,  a>  is  the  angular  frequency  in  rad/s,  9a  is  the  pitching  amplitude  and  (j)  is  the  phase  shift  angle 
of  the  pitching  motion  leading  the  plunging  motion.  The  reduced  frequency  k  ,  is  defined  as 


k  =  2nfc  t  u  =  coc  /  u  (3) 

where  /  is  the  flapping  frequency  in  Hertz  and  u  is  the  flow  velocity  of  the  far  field.  This  flapping  model  is 
illustrated  using  a  2D  airfoil  in  Fig.  3,  where  O  is  the  pivot  point  and  a  is  the  dimensionless  length  from  the  leading 
edge  to  the  pivot  point  with  respect  to  the  chord. 


Fig.  3  Illustration  of  the  flapping  model. 

The  Reynolds  number  (Re)  and  the  Strouhal  number  (St)  are 


Re =  puc / p 


(4) 


St  =  fA  lu  (5) 

where  p  is  the  fluid  density,  p  is  the  fluid  dynamic  viscosity  and  A  is  the  wake  width  and  can  be  estimated  using 
the  peak-to-peak  excursion  of  the  trailing  edge,  or  more  simply  by  twice  the  plunging  amplitude. 

In  classical  aerodynamics,  the  lift  coefficient  C,  ,  the  drag  coefficient  CD ,  and  the  pitch  moment  coefficient  CM 
can  be  defined  as 


CD  = 


0.5  pirS 


(6) 


Q  = 


0.5 puS 


(7) 


cM=- 


M  _ 


(8) 


0.5  pu2SL 

where  Fx  and  F  are  the  components  of  the  resulting  aerodynamics  force  along  horizontal  (parallel  with  u 
direction)  and  vertical  (normal  to  u  direction)  directions,  respectively,  S  is  the  reference  area  and  equals  to  cin 
value  for  a  2D  problem,  M _  is  the  pitch  moment  with  respect  to  the  pivot  point  and  L  is  the  reference  length  and 

equals  to  c  here.  The  time-average  thrust  coefficient  CT  and  the  power  input  coefficient  Cp  in  one  flapping  cycle 
can  be  calculated  by  Eqs.  (9)-(10).  Correspondingly,  the  propulsive  efficiency  is  defined  in  Eq.  (11). 


—  —  1  rt+T 

CT  =  ~CD  =  — — CD  (t)dt , 


(9) 


CP  =j{\"T  CL{t)Kt)dt  +  \"T  CM{t)0{t)cdt),  (10) 

7  =  CTu  /  Cp ,  (11) 

where  T  is  the  period  in  seconds  with  T  =  1  I  f  and  y(t )  and  9{t)  are  the  first  order  time  derivation  of  y(t)  and 
0(1) ,  respectively.  For  simplicity,  we  also  use  CD  to  stand  for  CD  ,  C,  for  C,  in  the  next  sections. 


B.  Two  Dimensional  CFD  Analysis 


The  commercial  CFD  solver  FLUENT  [17]  is  employed  to  simulate  the  unsteady  flow  fields  around  the  moving 
wings  with  predefined  motions.  The  time-dependent  Navier-Stokes  equations  are  solved  using  the  finite  volume 
method,  assuming  incompressible  laminar  flow.  The  mass  and  momentum  equations  are  solved  in  a  fixed  inertial 
reference  frame  incorporating  a  dynamic  mesh.  The  continuity  and  the  momentum  equations  are  given  as 

V-V  =  0 ,  (12) 

p\^-  +  V-Vv\  =  -Vp  +  pS/2V ,  (13) 

where  V  and  p  are  velocity  and  pressure,  respectively. 

The  hybrid  mesh  which  is  shown  schematically  in  Fig.  4  is  employed  to  compute  the  flow  field.  The 
computational  domain  is  divided  into  two  distinct  zones:  moving  zone  and  re-meshing  zone.  The  moving  zone 
consists  of  the  C-type  structured  quadrilateral  mesh,  and  the  re-meshing  zone  consists  of  the  unstructured  triangular 
mesh.  The  airfoil  is  located  in  the  center  of  the  computational  domain,  and  the  no-slip  wall  boundary  condition  is 
applied.  The  spatial  scale  of  each  zone  and  corresponding  boundary  condition  are  also  shown  in  Fig.  4.  The  whole 
moving  zone  mesh,  including  the  interfaces  between  these  two  zones,  moves  with  the  airfoil  together  according  to 
the  predefined  airfoil  motion.  This  means  re-meshing  only  occurs  at  a  distance  of  20  to  45  reference  lengths  away 
from  the  airfoil  body,  which  insures  that  the  flow  simulation  around  the  airfoil  is  somewhat  affected  by  the  moving 
mesh.  The  C-type  mesh  structure  in  the  very  close  neighborhood  of  the  airfoil  is  shown  in  Fig.  5,  and  the  grid  size  is 
201x101  nodes  (201  along  every  single  airfoil  surface,  101  in  the  vertical  direction)  with  the  thickness  of  the  first 
layer  grid  around  the  airfoil  equal  to  0.0002c.  The  hybrid  mesh  is  generated  in  GAMBIT  v2.3.16  [26]. 


Tn  pave  grid  n-rnnhing  aim r 


Fig.  4  Flybrid  mesh  topology  with  boundary  conditions. 


Fig.  5  C-type  grid  very  close  to  the  airfoil. 


To  validate  the  accuracy  of  the  present  approach,  simulations  are  performed  in  7  periods  with  500  and  1000  time 
steps  in  one  plunging  period  under  the  conditions  of  £=2.0,  h=0.4,  «=34.7  m/s  (Ma=0.1),  c=0.064  m,  and 
Re=\  ,0x  1 04  respectively.  The  coupling  between  the  pressure  and  the  velocity  is  achieved  by  means  of  the 
SIMPLEC  algorithm.  Meanwhile,  the  discretizations  of  pressure  and  momentum  terms  are  the  Second  Order 
scheme  and  the  Second  Order  Upwind  scheme.  The  time  discretization  is  the  First  Order  Implicit  scheme,  which  is 
a  more  straightforward  method  in  FLUENT  for  the  dynamic  mesh  module  [27].  The  time  variation  of  the  plunging 
position  and  the  time  histories  of  CD  in  one  period  are  shown  in  Fig.  6,  with  the  comparisons  with  results  obtained 
by  Tuncer  [28]  and  Miao  [22].  These  four  results  have  a  good  agreement,  though  different  mesh  schemes  are 
employed  in  these  studies.  Figure  6  also  shows  that  500  time  steps  in  one  cycle  are  good  enough  to  get  the 
concerned  details.  Based  on  the  validation,  201x101  size  grid  with  the  first  layer  thickness  of  0.0002c,  and  500  time 
steps  are  employed  for  all  of  the  simulations  in  the  next  sections. 


Fig.  6  Validation  for  the  plunging  motion. 


C.  Optimization  Formulation 


The  path  optimization  problem  for  a  NACA0012  airfoil  [7]  is  formulated  under  the  flight  conditions  of  k  =  1 , 
u  =  3.0  m/s,  c  =  0.05  m,  and  Re  =  le4  as  follows: 

Find  h,  60,</ S 

tominimize  f(h,  0o,</)  =  ( \~P)^T  ^max  +p  u  Vax  (14) 

^T,  min  ^T,  max  7  min  7  max 

subject  to  0.8  2  /i  <  2.5,  10.0  2  90  <  30.0,  80.0  <  tj>  <  135.0 

where  /?  is  the  weight  factor,  CT  min  and  CT  max  are  the  minimum  and  the  maximum  time  average  thrust  coefficients 
in  the  design  array  for  constructing  the  Kriging  model,  and  //mm  and  are  the  minimum  and  the  maximum 
propulsive  efficiency  in  the  design  array  for  constructing  the  Kriging  model.  This  objective  function  is  equivalent  to 
maximizing  the  time  average  thrust  coefficient  and  the  propulsive  efficiency  simultaneously  under  the  given  weight 
factors. 

D.  Surrogate  Model  with  the  Kriging  Method 

In  this  research,  the  Kriging  model  is  employed  as  the  approximation  method,  and  the  orthogonal  array  as  the 
design  of  experiments  technique.  The  Kriging  method  initiated  from  the  geostatistical  field  can  be  considered  as  a 
member  of  the  family  of  linear  least  square  estimation  algorithms  [11].  The  Kriging  model  combines  a  global 
polynomial  model  plus  a  localized  departure.  This  is  why  the  Kriging  model  is  expected  to  find  a  well  global 
approximation  at  an  unobserved  location  for  a  nonlinear  problem  [11,  12].  Based  on  the  design  space  described,  h  is 
divided  into  7  levels,  d0  5  levels,  and  ^  6  levels  as  shown  in  Table  1.  The  design  array  in  this  paper  consists  of  the 
union  of  the  49  samples  obtained  using  a  classical  orthogonal  array  and  the  8  vertices,  which  means  55  samples  total. 
The  high  fidelity  results  are  obtained  using  FLUENT  under  the  given  flight  conditions. 

In  order  to  make  the  surrogate  model  be  more  flexible,  two  Kriging  models  are  constructed  to  fit  the  time  average 
thrust  coefficient  and  propulsive  efficiency,  respectively.  Five  design  points  not  belonging  to  the  samples  are 
selected  to  evaluate  these  two  surrogate  models  as  shown  in  Table  2.  The  largest  relative  error  between  the  high 
fidelity  model  and  the  surrogate  model  is  5.3%,  and  these  two  Kriging  models  are  supposed  to  be  good  enough  to 
predict  the  time  average  thrust  coefficient  and  the  propulsive  efficiency. 


Table  1  Design  space,  factors  and  levels 


No. 

1 

2 

3 

4 

5 

6 

7 

h 

0.80 

1.00 

1.20 

1.50 

1.80 

2.00 

2.50 

0a  (deg) 

10.00 

15.00 

20.00  : 

25.00 

30.00 

<t  (deg) 

80.00 

90.00 

100.00  110.00 

120.00 

135.00 

Table  2  Validation  of  the  Kriging  models 

No. 

h 

00  (deg) 

fzi  (deg) 

cD 

(sModel) 

CD 

(FLUENT) 

Error 

t]  (sModel) 

//  (FLUENT) 

Error 

1 

1.60 

23.50 

103.40 

-1.399019 

-1.449604 

3.50% 

0.297454 

0.296856 

0.20% 

2 

1.36 

29.60 

97.80 

-1.159442 

-1.138866 

1.80% 

0.451900 

0.476958 

5.30% 

3 

1.73 

23.80 

100.70 

-1.554946 

-1.591267 

2.30% 

0.280245 

0.277590 

1.00% 

4 

1.52 

26.90 

87.20 

-1.201629 

-1.218762 

1.40% 

0.341126 

0.355931 

4.20% 

5 

1.55 

28.60 

94.90 

-1.423153 

-1.468366 

3.10% 

0.374416 

0.390918 

4.20% 

E.  Optimization  Results 

The  weight  factor  f> ,  is  set  to  be  0.5,  and  the  ga  function  embedded  in  MATLAB  [20]  is  employed  to  solve  the 
optimization  problem.  The  optimization  history  is  presented  in  Fig.  7.  Correspondingly,  the  time  average  thrust 
coefficient  and  the  propulsive  efficiency  are  1.10  and  0.4731,  and  the  optimum  design  point  is  h  =  1.3265  , 
Ga  =  30.0  ,  and  (j)  =  97.6624  .  It  should  be  mentioned  that  this  optimum  is  similar  to  Tuncer’s  optimum  with 


£  =  0.5  [7], 


Fig.  7  Convergent  history  of  the  genetic  algorithm. 

F.  Three  Dimensional  CFD  Analysis 

The  grid  near  the  3D  wing  is  shown  in  Fig.  8,  and  this  3D  hybrid  grid  is  validated  using  a  3D  wing  with  pure 
plunging  motion  under  the  conditions  of  £=3.64,  £=0.175,  u= 0.3  m/s,  c=0.1  m,  b(semi-span)=0.3  m,  and 
/?e=3.0x  1 04.  The  results  are  compared  with  Heathcote’s  experimental  results  [23]  and  Tang’s  computational  results 
[29]  illustrated  in  Fig.  9.  Both  of  the  present  results  and  Tang’s  results  have  a  little  derivation  from  the  experimental 
results.  Flowever,  the  present  results  match  well  with  the  previous  computational  results. 


Fig.  8  A  structured  grid  very  close  to  the  wing. 
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Fig.  9  Validation  for  a  wing  with  pure  plunging  motion. 


In  the  case  of  2D  computations,  the  wing  tip  vortexes  and  the  spanwise  flows  are  ignored  compared  to  the  3D 
simulations.  However,  we  still  can  expect  that  the  2D  results  are  an  acceptable  approximation  of  the  3D  results  [30]. 
A  wing  with  AR= 6  and  the  optimum  kinematic  parameters  are  simulated.  The  time  average  thrust  coefficient  and 
the  propulsive  efficiency  are  shown  in  Table  3.  It  is  observed  that  C7  and  t]  obtained  by  the  Kriging  model  and 
FLUENT  at  the  optimum  are  very  similar  (the  relative  error  is  less  than  5.5%).  This  validates  the  Kriging  model 
again.  The  optimum  of  the  2D  problem  is  assumed  to  be  the  optimum  kinematic  parameters  for  the  3D  wing, 
though  they  have  a  little  derivation  in  the  time  average  thrust  coefficient,  up  to  13.24%.  Topology  optimization  in 
the  next  sections  is  based  on  this  optimum. 


Table  3  Comparison  between  2D  and  3D  results 


No. 

h 

^0 

<t 

CV(  Kriging) 

CtCFLUENT) 

'/(Kriging) 

//(FLUENT) 

2D 

1.3265 

30.0000 

97.6624 

1.1000 

1.0705 

0.4731 

0.5001 

3D 

1.3265 

30.0000 

97.6624 

- 

0.9453 

- 

0.4992 

IV.  Topology  Optimization 


As  mentioned  earlier,  the  flapping  motion  is  described  by  the  pitching  motion  and  plunging  motion.  The  flapping 
motion  is  optimized  by  path  optimization.  The  unsteady  aerodynamics  analysis  is  performed  using  the  results  of 
path  optimization.  The  dynamic  pressures  of  the  surface  are  obtained.  The  pressures  are  applied  to  the  finite 
element  model  of  the  wing  for  topology  optimization.  Topology  optimization  is  used  to  find  the  optimal  lay-out  of  a 
structure  within  a  specified  region.  In  general,  the  objective  function  to  be  minimized  in  topology  optimization  is 
the  compliance  with  an  equality  constraint  upon  the  volume  fraction.  Therefore,  topology  optimization  is  performed 
using  the  results  of  path  optimization  in  order  to  reduce  the  weight  and  maximize  stiffness  to  find  the  distribution  of 
reinforcement. 

A.  Finite  Element  Model  of  the  Flapping  Wing 

The  surface  pressures  are  obtained  by  using  the  optimum  values  of  path  optimization.  First,  the  finite  element 
model  is  generated  based  on  the  grid  scheme  of  the  unsteady  aerodynamics  analysis.  However,  the  finite  element 
model  is  not  suitable  for  structure  analysis  because  the  aspect  ratio  of  the  elements  is  too  high.  Therefore,  the  grid  is 
changed  for  structural  analysis  and  Fig.  10  shows  the  finite  element  models.  All  the  element  pressures  of  Fig.  10  a) 
are  obtained  from  unsteady  aerodynamics  analysis.  The  pressures  are  mapped  to  the  finite  element  model  for 
structural  analysis  using  a  0-order  interpolation  as  illustrated  in  Fig.  10. 
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a)  b) 

Fig.  10  Application  of  the  pressures  of  the  finite  element  model:  a)  unsteady  aerodynamic  analysis  model,  b) 

structure  analysis  model. 


Figure  11  presents  the  finite  element  model  of  the  wing  for  topology  optimization.  In  the  flapping  wing,  the 
center  point  of  the  leading  edge  is  attached  to  the  frame  as  illustrated  in  Fig.  1.  The  boundary  condition  is  imposed 
on  the  13  nodes  of  the  center.  All  the  degrees  of  freedom  in  the  six  directions  are  fixed.  The  finite  element  model  is 
composed  of  two  shell  element  layers  as  illustrated  in  Fig.  12.  In  Fig.  12,  the  finite  element  model  is  divided  into 
the  design  domain  and  the  non-design  domain. 


Fig.  11  The  finite  element  model  of  the  wing  for  topology  optimization. 


Non-design  domain 


Design  domain 


Fig.  12  The  design  and  non-design  domain  for  topology  optimization. 


The  design  domain  is  the  inner  layer  of  the  shell  elements.  One  layer  is  composed  of  3,464  shell  elements. 
Common  nodes  are  used  for  the  outer  and  inner  shell  elements.  The  Young’s  modulus  of  the  design  domain  (carbon) 
is  170.3  GPa,  the  Poisson’s  ratio  is  0.3  and  material  density  is  2001.2  kg/m3.  The  Young’s  modulus  of  the  non¬ 
design  domain  (mylar)  is  3.28  GPa,  the  Poisson’s  ratio  is  0.4  and  material  density  is  1400.6  kg/m’.  The  thickness  of 
the  design  domain  is  1.0  mm  and  the  thickness  of  the  non-design  domain  is  0.  1  mm.  The  pressures  are  always 
applied  on  the  non-design  domain  and  the  non-design  domain  is  supported  by  the  design  domain.  Therefore,  the 
topology  of  the  reinforcement  (carbon)  is  changed  during  the  iteration  of  topology  optimization. 

In  order  to  perform  dynamic  topology  optimization,  the  surface  pressures  are  calculated  at  every  time  step  of 
dynamic  analysis.  The  pressures  are  applied  to  the  finite  element  model. 


B.  General  Formulation  of  Topology  Optimization 


Topology  optimization  determines  the  state  of  material  distribution  in  the  structural  domain.  The  intermediate 
variable  is  used  for  the  design  variables  in  the  density  method  and  has  a  value  from  0  to  1 .  These  design  variables 
determine  the  existence  of  the  materials.  An  element  does  not  exist  when  the  value  of  the  design  variable  is  zero 
and  exists  when  the  value  of  the  design  variable  is  1  [31].  The  general  purpose  of  topology  optimization  is  to 
maximize  stiffness.  The  objective  function  is  the  average  compliance  and  the  design  variables  are  the  intermediate 
density  variables.  The  topology  optimization  formulation  is  as  follows: 

Find  b 

to  minimize  compliance  (15) 


subject  to  Zveb<F 


where  the  design  variable  vector  b  is  the  density  of  elements,  ve  is  the  volume  of  each  element,  and  F  is  the  total 
volume  of  the  structure.  Mainly,  strain  energy  is  used  for  the  objective  function  as  compliance  in  linear  static 
topology  optimization  because  compliance  is  twice  as  large  as  strain  energy. 

The  design  variables  in  topology  optimization  do  not  directly  correspond  to  the  existence  of  the  material. 
Therefore,  the  material  property  update  method  is  utilized  to  use  the  results  of  the  previous  iteration  in  the  next 
iteration.  The  update  ride  of  the  material  property  of  each  element  is  as  follows  [31]: 


E  =  b<t>PEn 


(16) 


Pm=bw?Po  (17) 

E  is  the  Young’s  modulus  vector  for  all  the  elements  and  it  represents  the  stiffness  of  an  isotropic  elastic  material. 
p,„  is  the  vector  for  the  density  of  the  material.  The  value  of  E0  means  the  initial  Young’s  modulus  and  the  value 
of  p0  means  the  initial  density  of  the  material.  In  general,  the  value  of  p  is  2  or  3  and  the  value  of  q  is  1  [24].  This 
update  method  for  the  material  property  is  used  each  iteration  in  the  commercial  software  system.  Static  topology 


optimization  proceeds  this  way. 


C.  Formulation  of  Dynamic  Topology  Optimization 


All  the  forces  in  the  real  world  act  dynamically  on  structures.  The  surface  pressures  of  the  wing  are  changed  in 
accordance  with  the  time  during  the  flapping  motion.  Thus,  design  and  analysis  should  be  performed  based  on  the 
dynamic  loads.  Dynamic  topology  optimization  is  performed  to  find  the  distribution  of  reinforcement  on  the 
flapping  wing. 


Strain  energy  Strain  energy 


- Initial  model  - Optimized  model  .  Average  of  optimized  model 


a)  b) 

Fig.  13  The  graph  of  the  strain  energy:  a)  minimization  of  summations,  b)  minimization  of  the  peaks. 

In  dynamic  topology  optimization,  some  researchers  claimed  that  the  summation  of  the  strain  energy  is  minimized 
during  the  entire  time  [25,  32].  Dynamic  topology  optimization  is  carried  out  for  the  stiffest  structure  with  a  specific 
mass  in  the  time  domain.  Figure  13  shows  examples  of  the  strain  energy  between  the  initial  model  and  the 
optimized  model.  The  initial  model  has  a  100  percent  mass  of  the  structure  and  the  optimized  model  has  a  specific 
mass  according  to  the  constraint  on  the  mass  fraction.  Thus,  the  strain  energy  of  the  optimized  model  is  larger  than 
that  of  the  initial  model  under  the  same  dynamic  load. 

Figure  13  a)  shows  the  minimization  of  strain  energy  summation  during  the  entire  time  and  Fig.  13  b)  shows  the 
minimization  of  the  peak  values  of  the  strain  energy.  In  the  optimized  model,  the  strain  energy  of  Fig.  13  a) 
increases  rapidly  at  a  particular  time,  while  the  average  strain  energy  of  Fig  13  a)  is  smaller  than  that  of  Fig  13  b).  It 
means  that  the  deformation  of  the  structure  can  be  very  large  at  a  particular  time.  In  Fig.  13  b),  the  average  of  the 
strain  energy  is  larger,  while  the  peak  is  smaller.  This  means  that  the  deformation  of  the  case  in  Fig.  13  b)  is 
uniform  in  the  time  domain.  Since  a  large  deformation  for  even  a  short  time  is  dangerous  to  the  flapping  motion,  the 


peaks  of  the  strain  energy  are  minimized  in  dynamic  topology  optimization.  The  peaks  of  the  strain  energy  are 
found  by  dynamic  analysis.  As  the  design  proceeds,  the  peaks  of  the  current  iteration  may  not  be  the  peaks  in  the 
next  iteration.  Thus,  the  strain  energies  for  a  few  time  steps  near  the  peaks  are  minimized  as  illustrated  in  Fig.  13  b). 
During  the  included  time  steps,  the  summation  of  strain  energy  is  minimized.  The  peaks  of  strain  energy  are 
decreased  and  the  topology  of  the  structure  is  obtained  for  a  specific  mass. 

A  dynamic  optimization  method  has  been  developed  for  size  and  shape  optimizations  based  on  ESLSO.  ESLSO 
includes  dynamic  characteristics  without  complex  differential  equations.  The  efficiency  of  ESLSO  has  been  verified 
by  various  researches  [16,  33-37].  As  mentioned  earlier,  ESLSO  is  expanded  for  dynamic  topology  optimization  of 
the  flapping  wing.  Equivalent  static  loads  (ESLs)  are  made  to  generate  the  same  response  field  as  that  from 
dynamic  pressures  at  each  time  step  of  unsteady  aerodynamics  analysis.  At  the  time  steps  near  the  peaks  of  the 
strain  energy,  ESLs  are  calculated  and  ESLs  are  used  as  the  external  pressures  in  static  topology  optimization. 

Dynamic  topology  optimization  using  ESLs  is  formulated  as  follows: 

Find  bf,  /  =  0,  ...,3200 

to  minimize  peaks  of  strain  energy  (18) 

subject  to  mass  <  massinitial  x  0.2 

0<6,  min  <b,<\  .0 

The  design  variable  6,  is  the  design  variable  of  the  /th  element  in  the  design  domain.  The  objective  function  is 
the  summation  of  the  strain  energy  at  the  time  steps  near  the  peaks  of  the  strain  energy.  The  constraint  is  that  the 
mass  of  the  structure  should  be  less  than  20  percent  of  the  initial  mass. 

D.  ESLSO  for  Topology  Optimization 

ESLSO  is  a  method  for  dynamic  optimization  on  structures.  The  loads  that  are  made  to  generate  the  same 
displacement  field  as  dynamic  loads  at  each  time  steps  are  called  ESLs.  The  ESLs  consider  all  the  behaviors  of  a 
structure  through  time.  Then  the  ESLs  are  applied  to  static  optimization  because  the  loads  remove  the  time 
dependent  constraint  conditions.  Calculated  ESLs  at  each  time  step  are  used  as  multiple  loading  conditions  in 


topology  optimization.  The  process  of  calculating  ESLs  is  explained.  The  governing  equation  of  dynamic  analysis 


is 


M(b)z(t)  +  C(b)z(0  +  K(b)z(0  =  f(t);  t  =  (),■■■. I  (19) 

where  M(b)  is  the  mass  matrix,  C(b)  is  the  damping  matrix,  K(b)  is  the  stiffness  matrix,  z (t)  is  the  dynamic 

displacement  vector,  t  is  the  time  and  1  is  the  number  of  time  steps.  ESL  is  calculated  as 

feq(s)  =  K(b)z(  <,•)/  s  =  0,---,I  (20) 

where  f  (s)  is  the  ESLs  vector,  K(b)  is  the  stiffness  matrix  and  s  is  the  loading  case  which  is  equivalent  to  t  in 

Eq.  (19).  The  number  of  ESLs  is  the  same  as  that  of  the  time  steps. 

ESLSO  was  originally  developed  for  size/shape  optimization.  A  preliminary  study  on  topology  optimization  using 
ESLs  was  performed  for  a  small  scale  standard  problem  [25].  As  mentioned  earlier,  since  the  summation  of  the 
strain  energy  in  the  entire  time  domain  was  minimized,  the  optimum  topology  was  not  practical.  In  this  research, 
the  method  is  modified  to  minimize  the  peaks  and  applied  to  the  design  of  a  flapping  wing.  The  process  of  topology 
optimization  using  ESLSO  is  as  follows: 

Step  1.  Set  cycle  number  k  =  0  ,  and  the  design  variables  b(A)  =  b<0) . 

Step  2.  Perform  dynamic  analysis  with  b'7' 1 . 

Step  3.  Calculate  ESLs  in  the  time  domain  by  using  Eq.  (20). 

Step  4.  Perform  static  topology  optimization  as  follows: 

Find  b 

to  minimize  /( b,  z ) 


subject  to  K(b)z(s)  =fe?(^);  s=l,---,q 


(21) 


where  the  design  variables  are  the  density  of  each  element  and  the  objective  function  is  the  summation  of  strain 
energy  at  the  time  steps  near  the  peaks  of  the  strain  energy. 

Step  5.  When  k  =  0  ,  go  to  Step  6.  Otherwise,  check  the  convergence. 

If  the  convergence  criterion  is  satisfied,  the  optimization  process  terminates.  If  the  convergence  criterion  is  not 
satisfied,  go  to  Step  6. 

Step  6.  Update  the  design,  set  k  =  k  + 1 ,  and  go  to  Step  2. 

The  process  is  repeated  until  the  convergence  criterion  is  satisfied.  The  repetition  process  is  called  a  cycle.  The 
difference  of  dynamic  and  static  responses  is  decreased  as  the  repetition  proceeds.  The  overall  process  is  illustrated 
in  Fig.  14. 


Fig.  14  The  process  of  topology  optimization  using  ESLSO. 


1.  Convergence  Criterion 

A  convergence  criterion  can  be  defined  by  borrowing  the  existing  one  for  size  and  shape  optimizations.  It  is 
defined  as 


i|bw-b,lt_1,|<£1  (22) 

where  sl  is  a  user-defined  small  number.  Equation  (22)  is  checked  in  Step  5.  After  static  topology  optimization,  if 
the  change  of  design  variables  from  the  previous  cycle  is  small,  the  process  terminates.  However,  it  is  difficult  to 
satisfy  this  condition  in  topology  optimization  because  some  design  variables  have  a  repeated  value  of  0  and  1. 
Therefore,  a  new  convergence  criterion  is  defined  as 

countifi >  Sj)  <  «xs2  (23) 

where  s2  is  another  user-defined  small  number.  In  Eq.  (23),  the  number  of  design  variables  which  change  more 
than  6  is  counted  first.  If  the  counted  number  is  less  than  certain  value  nxs2,  then  the  process  terminates.  The 
convergence  criterion  in  Eq.  (23)  is  utilized. 

2.  Method  for  Design  Update 

In  general  static  topology  optimization,  the  update  method  of  the  material  property  is  in  Eqs.  ( 1 6)-(  1 7).  However, 
this  update  method  is  used  within  static  topology  optimization.  The  update  method  between  the  cycles  is  newly 
defined  for  dynamic  topology  optimization.  The  optimum  design  variables  from  static  topology  optimization  have 
the  values  in  the  continuous  apace  from  0.0  to  1.0.  After  static  topology  optimization  is  finished,  the  elements  less 
than  a  user  defined  small  number  c,  are  removed  and  the  design  variables  for  the  elements  larger  than  s3  are 
remained.  It  is  noted  that  s3  should  be  fairly  small  because  the  results  with  a  large  £\  deteriorate  the  connectivity 
of  the  elements.  Even  a  small  c,  removes  elements  quite  a  lot.  The  density  of  an  element  is  changed  so  that  a 
removed  element  has  0.0  and  the  remained  element  has  1.0.  The  new  density  is  updated  and  dynamic  analysis  is 
performed  with  the  new  densities  in  the  next  cycle. 


E.  Results  of  Topology  Optimization 


As  mentioned  earlier,  dynamic  topology  optimization  is  performed  to  find  the  reinforcement  of  the  flapping  wing. 
The  dynamic  pressure  from  the  3D  CFD  analysis  is  imposed  on  the  flapping  wing  as  the  external  loads.  The  peak  of 
the  strain  energy  is  minimized  in  dynamic  topology  optimization.  The  number  of  the  time  steps  near  a  peak  is  set  by 
five  and  ESLSO  is  used  for  dynamic  topology  optimization.  The  results  of  the  dynamic  topology  optimization  are 
illustrated  in  Fig.  15. 


a)  b) 

Fig.  15  The  results  of  dynamic  topology  optimization:  a)  Mass  fraction:  0.1,  b)  Mass  fraction:  0.05 

Figure  15  a)  is  the  result  of  dynamic  topology  optimization  of  the  flapping  wing  when  the  mass  fraction  is  0.1,  e, 
is  0.3,  E2is  0.2  and  s 3  is  0.00001.  Figure  15  b)  is  the  results  of  dynamic  topology  optimization  when  the  mass 
fraction  is  0.05,  8!  is  0.1,  s2  is  0.1  s}  is  0.00001.  The  remaining  elements  show  the  distribution  of  the 
reinforcement  of  the  wing. 

The  structure  of  the  flapping  wing  should  have  some  flexibility  as  well  as  certain  rigidity  and  a  small  amount  of 
rigidity  guarantees  flexibility  [38-40].  Flowever,  how  flexible  or  rigid  the  wing  should  be  is  not  fully  studied  yet  for 
each  flapping  wing.  This  research  only  shows  that  the  results  of  dynamic  topology  optimization  enable  to  keep  a 
certain  rigidity  with  a  specific  mass.  If  we  know  the  target  rigidity,  topology  optimization  can  be  performed  to  meet 
the  target  rigidity.  This  aspect  will  be  explained  in  the  next  section. 

V.  Future  Direction  for  the  Design  of  the  Flapping  Wing 

This  study  shows  the  possibility  of  optimizing  the  flapping  wing  system.  Flowever,  all  the  involved  disciplines  are 
separately  considered.  Therefore,  the  optimum  solution  is  obtained  in  a  local  sense.  In  the  future,  they  should  be 
considered  simultaneously.  A  direction  of  the  future  study  is  illustrated  in  Fig.  16.  To  include  the  stiffness  of  the 
structure,  path  optimization  should  be  simultaneously  performed  with  structural  optimization,  and  aeroelasticity 
analysis  is  required  in  this  process.  A  high  fidelity  structural  model  may  cost  too  much  in  this  optimization  process. 
A  low  fidelity  (simplified)  structural  model  can  be  exploited.  The  approximation  method  of  this  research  can  be 


utilized.  Then  the  target  for  rigidity  is  evaluated  and  dynamic  topology  optimization  is  performed  to  meet  the  target. 
The  external  loads  are  made  by  the  method  of  this  research.  For  a  final  detailed  design,  size/shape  optimization  is 
performed  to  meet  the  target  while  various  design  conditions  are  satisfied.  Eventually,  multi-disciplinary  design 
optimization  (MDO)  will  be  employed  because  the  disciplines  are  coupled. 


Path  and  structural  optimization  of  the  flapping  wing 

i 

k 

Aeroelasticity  analysis  with  a  simplified  structure 

Optimum  values  (path  and  target  rigidity) 

'  ’ _ 


Dynamic  topology  optimization  for  the  target 


t 


Size/  shape  optimization  for  the  structure 


Fig.  16  Scenario  for  future  direction. 


VI.  Conclusions 

A  flapping  system  is  optimized  and  path  optimization  is  conducted  to  enhance  the  performances  such  as  thrust  and 
propulsive  efficiency.  The  Kriging  method  is  employed  for  optimization  and  CFD  analysis  for  a  2D  airfoil  is 
utilized  for  the  optimization  process.  The  determined  path  is  utilized  for  dynamic  topology  optimization  of  the  wing 
structure.  The  peaks  of  the  strain  energy  are  minimized  to  find  the  stiffest  structure.  The  external  loads  are 
calculated  from  a  3D  CFD  airfoil  analysis.  Reinforcement  within  the  wing  is  determined  to  have  certain  rigidity 


with  minimum  mass.  Since  the  external  loads  are  dynamically  imposed  in  the  time  domain,  a  dynamic  topology 


optimization  method  is  developed  by  using  the  equivalent  static  loads.  Since  the  research  is  in  the  early  stage,  future 
direction  for  the  research  is  proposed. 
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